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Abstract 

I 

. 2 ^ wide variety of approaches to estimate the degree of synchrony between two or more spike trains have been proposed. One of the 
inost recent methods is the ISI-distance which extracts information from the interspike intervals (ISIs) by evaluating the ratio of 
^ the instantaneous firing rates. In contrast to most previously proposed measures it is parameter free and time-scale independent. 
^ However, it is not well suited to track changes in synchrony that are based on spike coincidences. Here we propose the SPIKE- 
■ ^ clistance, a complementary measure which is sensitive to spike coincidences but still shares the fundamental advantages of the 
^^ISI-distance. In particular, it is easy to visualize in a time-resolved manner and can be extended to a method that is also applicable 
to larger sets of spike trains. We show the merit of the SPIKE-distance using both simulated and real data. 
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1. Introduction 

\ One of the key challenges of neuroscience is to understand 
the neuronal code, the way information is represented in the 
neurons' spike trains. Determining whether different stim- 
uli can be distinguished from the spike train responses they 
trigger requires a notion of similarity of spike trains. This 
notion is typically expressed as a distance between spike 
trains. Reflecting the variety of neuronal coding hypothe- 
ses, several spike train distances with sensitivities to differ- 
ent features have been proposed (Victor and Purpura, 1997; 
Kreuz et al., 2007). Some of these distances follow the clas- 
sical view that firing rates play a central role in neural cod- 
ing (Barlow, 1972), and are thus based on the spike count 
(Victor and Purpura, 1996). However, in the last decades 
the fundamental relevance of temporal structure has been 
established (Theunissen and Miller, 1995), and accordingly 
many distances focus on the timing of spikes (Victor and 
Purpura, 1996; van Rossum, 2001; Schreiber et al., 2003). 
Special cases of temporal coding are addressed by methods 
that deal with, e.g., the detection of coincidences (Griin et 
al., 1999) and the identification of temporal spiking pat- 
terns (Abeles et al., 1988; Victor and Purpura, 1997). 

One of the most widely used spike train distances is 
the metric based on spike times introduced in Victor and 
Purpura (1996) which evaluates the cost needed to trans- 
form one spike train into the other using only certain el- 
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ementary steps. This method involves one parameter that 
sets the time-scale. In contrast, a more recent approach, 
the ISI-distance, is parameter free and time-scale adaptive 
(Kreuz et al., 2007). The ISI-distance is also easy to visu- 
alize in a time-resolved manner. On spike trains extracted 
from a simulated Hindemarsh-Rose network it could repro- 
duce a known clustering equally well as the best time-scale- 
optimized measure (Kreuz et al., 2007). Recently, two mul- 
tiple spike train extensions of the ISI-distance have been 
proposed, the local pairwise average and a multivariate 
measure based on the coefficient of variation (Kreuz et al., 
2009). 

The ISI-distance and its extensions are based on the 
relative length of simultaneous interspike intervals (ISI) 
and are thus well-designed to quantify similarities in the 
neurons' firing rate profiles. However, they are not opti- 
mally suited to track the type of synchrony which is me- 
diated by spike timing and in particular by changes in 
the fraction of coincident spikes. This specific kind of sen- 
sitivity is not only of theoretical importance but also of 
high practical relevance, since coincidences of spikes have 
been proven to be of high prevalence in many different 
neuronal circuits. Examples include olfactory bulb (Wil- 
son and Mainen, 2006), hippocampus (Best et al., 2001), 
somatosensory cortex (Romo and Salinas, 2003), auditory 
cortex (Schreiner et al., 2007), visual cortex (Usrey et al., 
1999; Priebe and Ferster, 2008), and retina (Meister and 
Berry, 1995; Field and Chichilnisky, 2007; Shlens et al, 
2008). 

Therefore, motivated by both the importance of tem- 



poral coding and the ubiquity of coincident neuronal fir- 
ing, we here propose the SPIKE-distance, a measure which 
uniquely combines the advantages of the ISI-distance - such 
as being parameter free, time-scale adaptive, and time- 
rcsolved - with a specific focus on spike timing. Given two 
spike trains, this focus is obtained by averaging for each 
time instant the absolute differences between the previous 
spike times and the following spike times and normalizing 
by the mean length of the interspike intervals. Furthermore, 
also for the SPIKE-distance we propose two extensions to 
the case of multiple spike trains, cither as a local pairwise 
average or as a multivariate measure. 

We use several sets of specifically designed and simu- 
lated spike trains in order to illustrate the properties of the 
SPIKE-distance and compare its performance to the one 
of previously published methods such as the ISI-distance, 
the Victor-Purpura distance and the correlation coefficient. 
In a first application to real data we employ the spiking 
activity of a population of retinal ganglion cells (Mcister 
and Berry, 1995; Field and Chichilnisky, 2007) as an ideal 
testing ground for the different measures. In particular, we 
evaluate their capability to reproduce the gradual decrease 
of synchrony between two retinal ganglion cells with the 
distance between their receptive fields (Meister et al., 1995; 
Shlens et al., 2006). As a second application to real data 
we employ the spike train distances to discriminate single- 
unit responses to taste stimuli recorded in the nucleus of 
the solitary tract (NTS) in rats (Di Lorenzo and Victor, 
2003; data available online at http://neurodatabase.org/, 
Gardner, 2004). In this case the performance of the various 
measures is quantified by the normalized mutual informa- 
tion of the confusion matrix (Abramson, 1963; Victor and 
Purpura, 1996). 

The remainder of the paper is organized as follows: In 
Section 2.1 we describe the new SPIKE-distance and its 
extensions. Subsequently, we use constructed examples to 
stress the different sensitivities of the ISI- and the SPIKE- 
distance (Section 2.2), and contrast their local behavior to 
the global properties of the Victor-Purpura distance (Sec- 
tion 2.3). We evaluate their capability to distinguish differ- 
ent levels of spike train synchrony (Section 3.1) and com- 
pare their performance on two real datasets (Section 3.2). 
Finally, we provide some details on the derivation of the 
SPIKE-distance A, and describe the previously introduced 
measures against which we compare (Appendix B) as well 
as the data that we use (Appendix C). In Appendix D we 
evaluate the performance of the new SPIKE-distance in the 
setups used in Kreuz et al. (2007, 2009). 

2. Methods 

In the following we will introduce the bivariate SPIKE- 
distance and its extensions. Like the ISI-distance (cf. Ap- 
pendix B.l) all these measures are based on instantaneous 
values, i.e., from the sequences of spike times we create 
time profiles with one value for each sampling point. The 



distances are then defined as the temporal average of the 
respective time profile. 

2.1. The bivariate SPIKE-distance and its extensions 

We denote with {t,!"'} = 4"^ •••i ^m'' ^^e spike times 
and with M„ the number of spikes for neuron n with n = 
1,...,N. 

2.1.1. Bivariate SPIKE-distance 

For each neuron n = 1, 2 we assign to each time instant 
the time of the previous spike 

4")(t)=max(t(")|t(") <t) i(")<t<i£, (1) 
and the time of the following spike 

4")(t)=min(t(")|t(")>t) 4"^<i<i2, (2) 

as well as the interspike interval 

x'i^lit)=t^^\t)-t'^\t). (3) 

We denote the instantaneous differences of previous and 
following spike times as 

Atp{t)=t'^\t)-t'^\t) (4) 

and 

AtF{t)=ti'\t)-t^^\t), (5) 

respectively. Note that by definition the intervals Atp{t) 
and Atp{t) do not cover the time instant t (see Fig. lA). 

An instantaneous spike time based measure of spike train 
distance is given by the regular mean of the two absolute 
differences of previous and following spike times. This func- 
tion is piecewise constant taking just one value within each 
interval of the pooled spike train. In order to be more lo- 
cal in time, we weight the two absolute differences \Atp{t)\ 
and |AtF(i)| such that the difference of the spikes that are 
closer dominate (Fig. IB). To this aim we denote with 

xP{t)=t-t^^\t) (6) 

and 

4")(t) = 4")(t)-t (7) 

the intervals to the previous and the following spikes for 
each neuron n = 1,2. The inverse values of their respective 
averages over both neurons (xp"''(t))„ and {x^\t))n serve 
as weights. 

Since the relevance of a certain spike time difference de- 
pends on the local rate (e.g., a small shift within a burst is of 
relatively higher importance than the same small shift be- 
tween two isolated spikes) . we normalize the average spike 
time difference by the mean interspike interval (here = 
2) 

1 ^ 

(41w)« = ^E41w- (8) 

n=l 

This way we also obtain time-scale invariance, i.e., the nor- 
malized spike time difference is the same for stretched or 
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Fig. 1. (color online) Illustration of the bivariate SPIKE-distance. A. 

in) 

Time instant t and its relation to the previous spikes tp and the 

following spikes tp"' as well as their differences. In this example the 
interspike interval of the first spike train is included in the interspike 
interval of the second spike train. The spikes of the two spike trains as 
well as their ISIs are shown in blue and green, respectively. The time 
instant under consideration is marked by a vertical dashed line, whereas 
horizontal solid arrows indicate the absolute differences |Atp(t)| and 
|Atp(t)| between the previous and the next spikes, respectively. B. 
Illustration of the weighting of the two contributions by the distance 
to the previous and the following spikes. Here the second spike of the 
second spike train comes first. The average distances to the previous and 
to the following spikes are represented by horizontal black lines. Note 
that for both examples we get the same value S'(t) = 0.367 (Eq. 9), 
since we just switched the following spikes which effected the individual 
intervals but not their sum. 

compressed spike trains. Using these quantities for locally 
weighted averaging and normalization we obtain the in- 
stantaneous spike time difference as (for a complete and 
more detailed derivation cf. Appendix A): 

|Atp(t)|(x^"^(t))„ + |AtF(i)l(4"^W)n 



S{t) 



(9) 



Finally, integrating over time yields the bivariate SPIKE- 
distance Ds: 

Ds^^ f S{t)dt. (10) 

The bivariate SPIKE-distance is bounded in the interval 
[0, 1]. The value is only obtained for perfectly identical 
spike trains. 

2.1.2. Averaged bivariate SPIKE-distance 

The averaged bivariate SPIKE-distance for a larger num- 
ber of spike trains is the bivariate SPIKE-distance averaged 
over all pairs of neurons. The same kind of time-resolved 



visualization as in the bivariate case is possible, because 
the two averages over time and over pairs of neurons com- 
mute. We thus can first calculate the instantaneous average 
S'^{t) over all pairwise instantaneous spike time differences 
^""(i) (Eq. 9) 
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and then average over time 



(11) 



S''{t)dt. 



(12) 



T jt=o 

Like the original SPIKE-distance, the averaged bivariate 
SPIKE-distance is restricted to the interval [0, 1]. 

2.1.3. Multivariate SPIKE-distance 

Calculating the average over pairs of spike trains can be- 
come very time-consuming as soon as the number of neu- 
rons N increases, since the computational cost scales with 
. In such cases a multivariate approach which scales only 
with N is computationally preferable. 

For each time instant we take the standard deviations 
over all previous and all following spike times and divide 
their mean by the mean instantaneous ISI. Similar to the 
bivariate case (cf. Section 2.1.1 and Appendix A), in order 
to be more local in time we weight the two terms such 
that the standard deviation of the spike times that are 
closer dominates. This way we obtain the instantaneous 
spike time standard deviation: 



S^'it) 



^[4")(i)]„(4"Ht))„ + a[4")(t)]„(4")(t)). 



(13) 

Finally, the multivariate SPIKE-distance Dg^ is obtained 
by integrating over time: 



Ds S^^{t)dt. 
^ Jt=o 



(14) 



The multivariate SPIKE-distance has as a lower bound 
but does not have an upper bound and in particular can 
attain values larger than 1. 

2.1.4. Practical considerations 

Since (xp"^(t))„ and {x'^\t))n depend explicitly on the 
time instant t, all three types of instantaneous SPIKE- 
measures S{t), S'^(i), and S"°(i) in Eqs. 9, 11, 13 are piece- 
wise linear rather than piecewise constant. Thus a new 
value has to be calculated for each sampling point and not 
just once per each interval in the pooled spike train. How- 
ever, this is only necessary when the localized visualization 
is desired. In case the distance value itself is sufficient, the 
short computation time can be even further decreased by 
representing each interval by the value of its center and 
weighting it by its length. This actually gives the correct 
result, the time-resolved calculation is a very good approx- 
imation for sufficiently small sampling intervals dt. 
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An alternative variant to the time- weighted average used 
in Eqs. 10, 12, and 14 is spike- weighted averaging where the 
SPIKE-distances are evaluated only after every new spike. 
This way each interval of the pooled spike trains contributes 
the same and not according to its length. For the sake of 
brevity in this study we restrict ourselves to presenting re- 
sults of the time- weighted average only, although the spike- 
weighted variant occasionally (for example for the data an- 
alyzed in Section 3.2.2) exhibits slightly better results. 

A technical detail concerns the ambiguous definition of 
the very first and the very last interspike interval as well 
as the initial distance to the previous spike and the final 
distance to the following spike. This issue can be resolved 
by placing for each spike train auxiliary leading spikes in 
the beginning of the recording at time t = and auxiliary 
trailing spikes at the end of the recording at time t — T. 
This leads to a systematic underestimation of the distance 
at the edges. While this limits the applicability of these 
methods to very sparse spike trains, the effect fades quickly 
with increasing numbers of spikes. A similar procedure has 
been used in the metric £)intcr™i:fix Victor and Purpura 
(1997). 

More details on the implementation as well as the 
Matlab source code for calculating and visualizing both 
the SPIKE- and the ISI-distances can be found under 
www.fi.isc.cnr.it / users / thomas.kreuz / sourcecode.html. 

2.2. Motivation: Sensitivity to spike coincidences and the 
spike train of origin 

The main difference between the SPIKE-based meth- 
ods and the previously proposed ISI-based methods (cf. 
Appendix B.l) is their different sensitivity to spike coin- 
cidences. The multiple spike trains shown in Fig. 2 were 
specifically constructed to illustrate this. In the first half 
(< 400 ms) of Fig. 2a all 20 spike trains are regular but with 
constant phase lags, while in the second half (400 — 800 ms) 
all spike trains are identical. Neither the averaged bivari- 
ate ISI-distance Df nor the multivariate ISI-distance Df^ 
distinguish between the intervals 100 — 300 ms (the first 
interval without the transients) and 400 — 800 ms, since 
the instantaneous frequencies are the same thus yielding a 
value of zero in both cases. The SPIKE-distances and 
on the other hand do distinguish, and higher values are 
obtained for the first interval in which spikes are maximally 
dispersed, whereas the coincident spikes in the second in- 
terval yield zero values indicating perfect synchrony. 

In the example of Fig. 2b we start with spike trains that 
are uncorrelated and basically random. At 200 ms and at 
500 ms two events are inserted, the first one perfectly syn- 
chronous and the second one with a certain amount of jit- 
ter. Subsequently, the spike trains become more and more 
regular which is achieved by creating events every 100 ms 
with decreasing amounts of jitter and closing with a per- 
fectly synchronous spiking event at 1900 ms. 

The perfectly synchronous spiking event embedded 
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Fig. 2. (color online) Instantaneous ISI- and SPIKE-measures for con- 
structed spike trains. On top we depict the spike trains and the mean 
instantaneous ISIs according to Eq. 3. Below that we show the four 
instantaneous measures (from top to bottom): the averaged bivariate 
ISI-ratio I^{t) (Eq. B.3), the multivariate coefficient of variation of the 
interspike intervals /'"(t) (Eq. B.5), the spike time difference S"'(t) (Eq. 
9), and finally the spike time standard deviation 5™(t) (Eq. 13). For 
the latter dashed vertical lines mark the occurrence of perfectly syn- 
chronous events. A. In the first interval all spike trains are regular, but 
with constant phase lags, then all spikes are regular with zero phase 
lag, i.e., identical. B. Spiking events with increasing jitter plus one noisy 
and one perfectly synchronous spiking event towards the end. For the 
instantaneous SPIKE-measures we also show the moving averages (of 
order 50 sampling intervals) at the bottom of this subplot. 

within the noise (200 ms) is recognized only by the two 
instantaneous SPIKE-measures S'^ and S"^ which decrease 
to zero. On the other hand, the instantaneous ISI-measures 
even show peaks which are caused by the normalization. 
Similar peaks occur at the second unreliable event (500 
ms). The SPIKE-measures on the other hand mark this 
event with a Mexican hat shaped event kernel ^ . The 
reason for the Mexican hat shape of these event kernels 
is intuitive: Right before the spiking event all differences 
between the following spike times are small as are all dif- 
ferences between the previous spikes right after the event. 
This affects respectively the first and the second half of 



^ Here the expression " Mexican hat" is not meant to strictly describe 
the difference between two Gaussians, but is rather used in a general 
sense to describe the succession "plateau - small drop - high peak - 
small drop -plateau" . 
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the numerator entering the instantaneous dissimilarities 
(of. Eq. 9). Accordingly, values are decreased with respect 
to the mean level. Furthermore, within the spiking event 
there is a small interval for which some of the spikes are still 
following spikes while others are already previous spikes. 
Thus for this small interval the respective differences be- 
tween previous and following spike times are very large 
and so is the value of the instantaneous SPIKE-measures. 

Subsequently, when the noise is gradually decreased, all 
measures, starting from a plateau with moderate fluctua- 
tions reflecting the initial inherent randomness, exhibit a 
consistent decrease. However, this drop occurs in a different 
manner for the instantaneous 1ST and SPIKE-measures. 
While the ISTmeasures decrease rather gradually with only 
a few smaller elevations, the SPIKE-measures mark the 
events with a series of peaks. These peaks are no longer 
Mexican hat shaped because between the events there is 
no noisy background level to decline to. With the increas- 
ing reliability of the spike events the peaks are getting 
more and more prominent. Consequently, a perfectly syn- 
chronous event like the one at 1900 ms could be represented 
as a peak with infinitely small width and maximum am- 
plitude. For such events the SPIKE-measures S'^ and S*™ 
yield the value zero (as expected for perfectly synchronous 
spike trains). The same value is also obtained for the in- 
tervals between two perfectly synchronous events. In order 
to distinguish these two cases we mark the occurrences of 
the perfect events by adding vertical dashed lines to the 
temporal profiles (see second half of Fig. 2A). 

The instantaneous SPIKE-measures are sensitive to 
spike coincidences and assign specific signatures to their 
time of occurrence. However, in case the focus of attention 
lies on the long-term changes in spike train synchrony, an 
appropriate moving average eliminates these short time- 
scale signatures and a gradual decrease similar to the one 
of the ISI-measures is obtained (cf. the two lowest subplots 
of Fig. 2B). 

The SPIKE-distances are not only sensitive to the timing 
of the individual spikes, but also to their spike train of 
origin (Fig. 3). To show this, we start with one spike train 
consisting of 10 equally distributed bursts of 10 spikes each. 
From this spike train we construct two sets of TV = 10 
spike trains with Af„ — 10 spikes each. This is achieved by 
assigning the individual spikes to their spike train of origin 
which is the inverse operation to pooling spikes from several 
spike trains. In the first set the spikes are distributed such 
that we have 10 unreliable events with one spike in each 
spike train (Fig. 3A), while in the second set the spikes are 
assigned to random spike trains keeping only the number 
of Mn = 10 spikes per spike train constant (Fig. 3B). As 
can be seen in the four lower traces, both the 1ST and the 
SPIKE-distances can clearly distinguish these qualitatively 
different behaviors. This is in contrast to the Peri-Stimulus 
Time Histogram (PSTH) which, like all other measures of 
spike train synchrony based on the PSTH or on the pooled 
spike train in general, is invariant to shuffling spikes among 
the spike trains. 
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Fig. 3. Two sets of spike trains demonstrating the sensitivity of the 
instantaneous ISI- and SPIKE-measures to the spike train of origin. By 
construction the pooled spike train of the two sets is identical consisting 
of 10 regularly distributed bursts. Only the distribution of the spikes 
among the individual spike trains differs. A. High reliability: Spikes from 
each burst are equally distributed among the spike trains. B. Irregular 
behaviour: All spikes are randomly assigned to the individual spike trains. 
Whereas the PSTH is by construction exactly identical in both cases, 
the ISI- and the SPIKE-distance can distinguish these cases resulting in 
low distances in the first case and in higher distances in the second case. 

2.3. Properties: Local behavior 

In order to illustrate the local behavior of the time- 
resolved ISI- and SPIKE-distances and compare it with 
the global behavior of the Victor-Purpura distance (cf. Ap- 
pendix B.2), we start with a simple example (Fig. 4) com- 
posed of two identical spike trains with just three spikes 
each. Then, while all other spikes are kept fixed, the inner 
spike of the second spike train is shifted in time. 

Exemplary spike trains and the results for the time- 
resolved bivariate measures I{t) and S{t) are shown in Fig. 
4A. Both measures show local signatures that are somehow 
counter-intuitive: They assign larger positive deviations for 
more closer spikes. The reason is as follows: As the dif- 
ference between the two inner spikes increases, the value 
of the instantaneous dissimilarities in the interval between 
these spikes (the peak height of the event kernel) decreases, 
however, at the same time the support interval of the peak 
value is getting longer as are the dissimilarity values in the 
remaining intervals (the drop height of the event kernel) 
although their support is getting smaller. Initially the sum 
effect is positive (compare the dark and the bright areas), 
however, as the difference between the two spikes is getting 
even larger, the relative importance of these two opposing 
effects changes. Accordingly, the corresponding distances 
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Fig. 4. (color online) Local behaviour of SPIKE-distances. A. All spikes 
are set to fixed positions except for the inner spike of the second spike 
train which is shifted relative to the inner spike of the first spike train. 
In this example this spike is shifted from 12 ms (solid lines) to 14 ms 
(dashed lines). Following the spike trains and the mean instantaneous 
ISI-values we depict the instantaneous ISI-ratio I{t) and the instanta- 
neous spike time difference S{t). The positive and negative deviations 
between these two profiles are marked by bright and dark shaded areas, 
respectively. In this example these two opposing contributions sum up 
to an increase of the distances. B. Dependence of the IS! distance Dj, 
the SPIKE-distance Dg and the Victor-Purpura distance Dy (for cost 
parameters cy equal to 0.01 and 0.2) on the relative shift of the two 
inner spikes. The relative positions of the spikes in A are marked by 
vertical lines. 

Di and Z?s shown in Fig. 4B exhibit an initial increase 
which then leads to a maximum value for an intermediate 
shift followed again by a decrease. As the different slopes 
of the initial increase in distance for positive and negative 
shifts prove, the relative importance of the two opposing 
effects does not only depend on the shift between the two 
intermediate spikes, but is also influenced by the relative 
lengths of the outer intervals. 

One peculiarity can be observed for the ISI-distance Z?i, 
a local minimum for a shift of 20 ms. This is the point 
for which the two spike trains are symmetric under time 
reversal. Accordingly, in the central interval between the 
two inner spikes the instantaneous measure I{t) equals zero 
and this leads to the low value of Di. 

In Fig. 4B we also show results of the Victor-Purpura 
distance for two cost values cy — 0.01 and cy = 0.2. 
A linear increase can be observed with the slope of the 
increase given by the cost parameter. Note that no distance 



values higher than 2 can be obtained, since for high costs 
and large differences of the two inner spikes shifting gets 
too expensive and it becomes more affordable to delete and 
to insert the spike for a cost of 2. 

The difference between the instantaneous measures and 
the Victor-Purpura distance stems from the fact that the 
first are local and the latter is global. For the SPIKE- 
distance the inner spike of the second spike train is ap- 
proaching the third spike of the first spike train which leads 
to the effect that the two spike trains are getting more sim- 
ilar again. The Victor-Purpura distance on the other hand 
matches the two third spikes and thus considers only the 
difference between the two inner spikes although they are 
further apart than the later spikes. 

3. Results 

3.1. Simulation: Capability to distinguish different levels 
of spike train synchrony 

As a first quantitative test for the different measures 
we evaluated whether they are able to track continuous 
changes in synchrony when the transition to synchrony is 
based on coincidences of spike times. Measures comprise 
the bivariate and multivariate 1ST and SPIKE-distances 
Df, D™, Dg, and D™, as well as the Victor-Purpura dis- 
tance Dy and the correlation coefficient C (cf. Appendix 
B.3). The latter two measures are both represented by two 
parameter values, one of which was optimized over differ- 
ent values of the respective parameter. 

We generated 25 spike trains with 100 spikes each. The 
spikes were randomly distributed within the interval [0, 1] 
using a sampling interval of 0.0001. In order to create spike 
trains with increasing levels of correlations we followed a 
correlation scheme introduced in Kreuz et al. (2006) and 
defined a matching parameter m that governed the fraction 
of shared spikes for each pair of neurons: for matching m = 
all spikes are randomly distributed, whereas for matching 
ni = 1 all spike trains are identical. 

Average values and standard deviations over 100 realiza- 
tions are shown in Fig. 5A. All distances equal for a perfect 
matching (to = 1). For decreasing matchings all average 
values increase monotonically, however, the shape of the in- 
crease differs. The multivariate methods first show a very 
steep and then a rather moderate increase. The averaged 
bivariate methods exhibit a more gradual increase except 
for the averaged bivariate ISI-distance whose increase is 
rather slow for low matchings. Furthermore, as indicated by 
the standard deviations of the distributions, the averaged 
bivariate methods are characterized by considerably nar- 
rower distributions. The average bivariate SPIKE-distance 
Z?g, the optimized Victor-Purpura distance IJy, and the 
optimized correlation coefficient C are the only measures to 
combine a rather constant increase with a low variability. 

Similar to the analysis of Kreuz et al. (2009) (see also 
Appendix D.2) we quantified the capability of the measures 
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Fig. 5. (color online) Capability to distinguish different levels of spike 
train synchrony. A. ISI- and SPIKE-distances as well as Victor-Purpura 
distance (normalized to its maximum value) and correlation coefficient 
versus the matching parameter. Values shown are averages and stan- 
dard deviations of distributions from 100 realizations with N = 25 spike 
trains and AI = 100 spikes each. For clarity we only show the stan- 
dard deviations for matching m = 0, the other standard deviations de- 
crease rather monotonously with the matching. B. Kolmogorov-Smirnov 
statistics showing the distinguishability of distributions with neighbor- 
ing matchings. Values for the different measures are plotted against the 
mean matching of the respective two distributions. C. Comparison of 
measures: Kolmogorov-Smirnov statistics averaged over all neighboring 
matchings. The vertical line separates the averaged bivariate and the 
multivariate measures. 



to distinguish distributions with neighboring matchings us- 
ing Kolmogorov-Smirnov statistics (Fig. 5B). For the aver- 
aged bivariate SPIKE-distance as weU as for the optimized 
optimized Victor-Purpura distance and the optimized cor- 
relation coefBcient the maximum value of 1 is obtained 
indicating pairs of non-overlapping neighboring distribu- 
tions. For the other measures deviations from 1 indicate 
overlapping distributions. This is quantified by the aver- 
age of the Kolmogorov-Smirnov statistics over all matching 
transitions (Fig. 5C). Perfect scores are obtained for the 
averaged bivariate SPIKE-distance, the optimized Victor- 
Purpura distance, and the optimized correlation coefficient 
followed by the averaged bivariate ISI-distance. In general, 
the multivariate measures yield lower averages with the ISI- 
distance again lagging behind the SPIKE-distance. 



3.2. Application to real data 

In Section 3.2.1 we investigate the spike train similar- 
ity within a complete population of retinal ganglion cells, 
while in Section 3.2.2 we evaluate how well the various 
measures can discriminate single-unit responses to differ- 
ent taste stimuli. Since in both applications we use se- 
tups which involve only similarities between pairs of spike 
trains, we only compare the basic bivariate measures and 
not the averaged bivariate or multivariate extensions. Thus 
the measures are narrowed down to the ISI-distance Di, the 
SPIKE-distance Ds, the Victor-Purpura distance Dy ^ , 
and the correlation coefficient C. In both applications for 
Dy and C we sampled the relevant range of the respective 
time-scale parameter equidistantly on a logarithmic scale. 

3.2.1. Spike train similarity among retinal ganglion cells 

We first applied the bivariate spike train distances to 
responses from an almost complete population {N — 105 
neurons) of ON parasol retinal ganglion cells (Shlens et al., 
2006, 2009; for more details on the data confer Appendix 
C.l). These recordings were performed using both white 
noise stimulation and constant, spatially uniform illumina- 
tion. The white noise stimulation was also used to identify 
the receptive field of the retinal ganglion cells by means 
of reverse correlation Chichilnisky (2001). As can be seen 
in Fig. 6A, the receptive fields of the complete population 
form an orderly mosaic that tiles visual space with minimal 
overlap. An exemplary snapshot of the spiking activity of 
the population responding to a white noise stimulation is 
shown in Fig. 6B. 

From this mosaic we extracted the pairwise distances be- 
tween the centers of the respective fields. The amount of 
synchronized firing between pairs of RGCs of the same type 
declines systematically with distance between the two cells 
(Mastronarde, 1983; Meister et al., 1995; De Vries, 1999; 
Shlens et al., 2006; Greschner, 2010). This weU described 
behavior was used as a benchmark against which we tested 
the capability of the various spike train distances to assign 
higher values to more distant cell pairs (or lower values to 
more distant cell pairs in case of the correlation coefficient 
C). Note that the distance between cell pairs does not ex- 
plain all the variance in the amount of synchronized firing 
and that this behavior was mainly studied using measures 
related to the correlation coefficient. 

For each of these measures and for both kinds of stimu- 
lations we first calculated all pairwise spike train distances 
for the complete population. Exemplary scatter plots of the 
measures versus the distance between the receptive fields 
are shown for one segment of the white noise stimulation 
in Fig. 7. All measures exhibit a pronounced saturation for 
higher distances between RFs but show a very clear depen- 



^ Note that among the time-scale dependent spike distances we 
restrict ourselves to the Victor-Purpura distance, since this was the 
measure that proved to be the best performer in previous measure 
comparisons (Kreuz et al. (2007, 2009), see also Appendix D). 
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Fig. 6. A. Receptive fields (Gaussian fit outlines) of an almost complete 
population of N = 105 ON-Parasol retinal ganglion cells. This is a 
reanalyzed subset of the first dataset of Shiens et al. (2006, 2009). B. 
Spike trains from these 105 cells during the first two seconds of one 
realization of white noise stimulation. 
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Fig. 7. Scatter plots of the ISI-distance Di, the SPIKE-distance Ds, 
the Victor-Purpura distance Dy (for logcy = 15), and the correlation 
coefficient C (for logb = —2) versus the distance between the centers 
of the receptive fields (RFs). Note that C is a measure of similarity 
and thus exhibits a decrease with the distance between the RFs. For 
both Dy and C we show the plots for the parameter value that yielded 
the highest absolute Spearman's rank correlation coefficient R (cf. Fig. 
8). The values of R for this example were 0.23 (Di), 0.36 (Ds), 0.17 
(Dv), and 0.39 (C), respectively. 

dency for lower distances. However, the measures differ in 
how pronounced this dependency is relative to the overall 
variance. While it is rather low for the Victor-Purpura dis- 
tance and the ISI-distance, it is much larger for the SPIKE- 
distance and the correlation coefficient. 

For all measures we quantified the relation between the 
spike train distances and the distances between the recep- 
tive fields using the absolute Spearman's rank correlation 



Fig. 8. Comparison of the absolute Spearman's rank correlation coeffi- 
cients (average and standard deviation) over 18 segments (10 s each) 
for both white noise stimulation and constant, spatially uniform illumi- 
nation. Measures comprise the ISI- and the SPIKE-distance, the Vic- 
tor-Purpura distance, and the correlation coefficient. The parameters of 
the latter two measures are varied on a logarithmic scale. 

coefficient R (Fig. 8). For both the white noise stimulation 
and the constant, spatially uniform illumination the high- 
est average i?-values are obtained for the optimized corre- 
lation coefficient C. This is closely followed by the SPIKE- 
distance Dg which in turn yields a considerable improve- 
ment over the ISI-distance Di. For the Victor-Purpura dis- 
tances Dy despite the optimization the highest coefficient 
is considerably lower than the values obtained for the other 
methods. Also results are not very robust since a more de- 
tailed investigation revealed that for different segments the 
actual maximum i?-value was attained for cost values cov- 
ering two entire logarithmic decades. 

3.2.2. Taste stimuli discrimination from single-unit 
recordings 

As a second application to real data we analyze one com- 
plete dataset from "neurodatabase.org" (Gardner, 2004). 
It consists of single-unit responses of three neurons (#4, 
^9, and #11) in the nucleus of the solitary tract (NTS) of 
rats to four different taste stimuli (Di Lorenzo and Victor, 
2003). For the three different neurons we have Nji = 19, 
23, and 16 repetitions, respectively, for each stimulus. For 
more details on the data confer Appendix C.2. 

For all measures and for each neuron separately we calcu- 
late the pairwise distance matrix for all 4 • realizations. 
Assuming stimulus-driven synchronization, a good discrim- 
ination should yield low distance values for responses to 
the same stimuli and high distance values to responses for 
different stimuli (and vice versa for the correlation coeffi- 
cient). 

In order to evaluate the separation of the spike trains 
into 4 different clusters, we follow Victor and Purpura and 
define the distance between the spike train Si and the clus- 
ter Cq, as {d{Si, Sj))a, where {■)a denotes the average over 
all spike trains in the cluster Cq.. From these distances we 
compute the normalized confusion matrix (Abramson, 
1963; Victor and Purpura, 1996) whose entry pc,^ is defined 
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Fig. 9. Normalized confusion mutual information for three neurons re- 
sponding to different taste stimuli. The same measures as in Fig. 8 were 
evaluated, only the range of the time-scale parameters was adapted in 
order to cover the maxima of discrimination for this task. 

as the probability that Cp is the closest cluster to a spike 
train belonging to Ca- For a perfect clustering the confu- 
sion matrix is diagonal, whereas each misclassification in- 
crements a non-diagonal element. The performance of the 
measures in discriminating the different stimuli is finally 
quantified by the normalized mutual information 



-fc = X! P-^P 



where Pa = J2bPab, and Pp = Y.bPbl3 and Imax denotes 
the maximum mutual information obtained for a correct 
classification. 

Results for the three neurons recorded during the taste 
stimuli discrimination task (Di Lorenzo and Victor, 2003) 
are shown in Fig. 9. For all three neurons the highest dis- 
crimination is obtained for the Victor-Purpura distance, 
however, as already reported in Di Lorenzo and Victor 
(2003) for different neurons the best discrimination is ob- 
tained for different coding schemes. The neurons ^4 and 
#11 exhibit their maximum for positive values of the cost 
parameter indicating that a temporal coding scheme allow 
for greater discrimination of the four stimuli. Only for neu- 
ron #:9 the curve does not increase considerably for pos- 
itive cost parameters. In this case the spike timing is less 
relevant and a rate coding scheme is sufficient for a proper 
discrimination of the tastants. Slightly below the maximum 
/c-values obtained for the Victor-Purpura distance follow 
the time-scale adaptive ISLdistance and SPIKE-distance. 
In this task there seems to be no clear distinction between 
their discrimination values regardless of the neuron's cod- 
ing scheme. Finally, for all three neurons the lowest dis- 
crimination capabilities are obtained by the correlation co- 
efficient. 

4. Discussion 

Motivated by the prevalence of coincident spiking in 
many different areas of the brain, we propose a SPIKE- 
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distance which, together with its multi-neuron extensions, 
can be applied to quantify the (dis) similarity of two or 
more spike trains. These measures are sensitive to the 
timing of spikes and are thus complementary to the ISI- 
distance (Kreuz et al., 2007) and its extensions (Kreuz et 
al., 2009) which are instead sensitive to the relative lengths 
of simultaneous interspike intervals. 

At the same time the SPIKE-distance and its extensions 
share many distinct properties with the ISI-distance, in par- 
ticular, they are parameter free and time-scale adaptive. 
The latter property can be preferable in applications to real 
data for which there is no validated knowledge about the 
relevant time scales. While time-scale dependent measures 
yield a functional characterization of the spike trains, a sin- 
gle valued method gives a more objective and comparable 
estimate of neuronal variability. Other drawbacks of time- 
scale dependent measures are the computational cost and 
the time and effort that is needed to find the right param- 
eter. Moreover, it is not at all guaranteed that there exists 
a right parameter. For example, spike trains that include 
different time-scales such as regular spiking and bursting 
might result in misleading conclusions, since any fixed pa- 
rameter will misrepresent either one of those dynamics. 

One of the main arguments for the use of a time-scale 
dependent measures is their potential insight into the pre- 
cision of the neuronal code (Victor and Purpura, 1996). 
A very common setup in experimental neuroscience is the 
recording of single- or multi-unit spike train responses for 
the repeated presentation of different classes of stimuli (for 
an overview see Victor and Purpura, 2005). The pairwise 
distance matrices over all realizations are calculated and 
used as input for a cluster analysis of the type used in Sec- 
tion 3.2.2. The time-scale for which the responses to dif- 
ferent stimuli are best discriminated is then assumed to be 
the precision of the neural code. However, a recent investi- 
gation of this kind of precision analysis (Chicharro et al., 
2011) has shown that the optimal time-scale obtained from 
the cluster analysis is far from being conclusive, but rather 
results in a non-trivial way from the interplay of different 
factors. These factors include the distribution of the infor- 
mation contained in different parts of the response and the 
degree of redundancy between them. For dynamic stimuli 
the optimal time-scale also depends on the stimuli statis- 
tics, in particular, the temporal distribution of the stimulus 
features to which the neurons respond. 

Despite the problems in the interpretation of the opti- 
mal time-scale, the Victor-Purpura distance is designed to 
test for specific codes ranging from a rate code to a coin- 
cidence detector (Victor and Purpura, 1997). The corre- 
lation coefficient focusses purely on the coincidences and 
thus reflects the most specific coding assumption. On the 
other extreme, the most general approach is taken by the 
ISI-distance and the SPIKE-distance which, in contrast to 
the Victor-Purpura distance, are time-scale adaptive and 
parameter-free. Thus the different measures can be ordered 
from measures tied to a specific hypothesis about the code 
(correlation coefficient) via intermediate measures (Victor- 



Purpura distance) to measures that are very general (isl- 
and SPIKE distance). Measures on different ends of this 
scale are complementary in nature. The ones capable of 
quantifying the degree of dissimilarity of the spike trains 
without relying on specific coding hypotheses are very well 
suited for an exploratory analysis. Oppositely, once a par- 
ticular coding scheme is hypothesized more specific mea- 
sures are needed for a confirmatory analysis. 

In the application to the retina data presented in Sec- 
tion 3.2.1 the correlation coefficient, which is focusscd on 
spike coincidences, yields best results, i.e., the highest abso- 
lute Spearman's rank correlation coefficients. On the other 
hand, in Section 3.2.2 the same measure fails to discrimi- 
nate the different taste stimuli. The former result is consis- 
tent with the prevalence of synchronized firing in the retina 
(Shlens et al., 2008), while in the latter taste discrimina- 
tion task the distinction seems to rely on coding schemes 
other than coincidences (Di Lorenzo and Victor, 2003). The 
more general ISI- and SPIKE-distances show a very adap- 
tive behavior and yield a good performance in both appli- 
cations. It should be noted that while both measures dis- 
criminate the taste stimuli about equally well, the SPIKE- 
distance outperforms the ISI-distance on the retina data. 
This is consistent with the previous observation: When syn- 
chronized firing is prevalent, the sensitivity to spike coinci- 
dences exhibited by the SPIKE-distance becomes relevant 
and leads to a considerable improvement with respect to 
the ISI-distance. 

Both ISI- and SPIKE-distance are conceptually simple, 
computationally efficient and easy to visualize in a time- 
resolved manner. This kind of visualization can be achieved 
because the distances are defined as temporal averages over 
time profiles, i.e., they are based on instantaneous values 
which are calculated from the sequences of spike times. 
However, it should be noted that each instantaneous value 
cannot be interpreted without looking at its local context. 
The SPIKE-distances are time-resolved only in the sense 
that they assign to each spiking event an event kernel. In 
case the focus of attention is on the long-term behavior, 
these local event kernels can easily be eliminated by ap- 
propriate moving averaging (as shown in Fig. 2b). Another 
caveat regarding all measures of spike train synchrony is 
that they generally assume a zero time delay between spike 
trains. Accordingly, in case of a non-zero time delay, e.g., 
caused by finite signal transmission times between different 
brain regions, it should be detected and eliminated before- 
hand (Waddefi et al, 2007; Nawrot, 2003). 

The measures used in this study quantify the level of 
synchrony within one set of two or more spike trains. For 
the Victor-Purpura and the van Rossum distances (van 
Rossum, 2001, cf. Appendix D) there also exist extensions 
that estimate the synchrony between two populations of 
neurons (Aronov, 2003 and Houghton and Sen, 2008, re- 
spectively) . Corresponding extensions for both the ISI- and 
the SPIKE-distance will be presented in a forthcoming 
study (Kreuz et al., 2011). 

We close by pointing out once more that the Mat- 



lab source code for calculating and visualizing both the 
1ST and the SPIKE-distances can be downloaded at 
www.fi.isc.cnr.it/users/thomas.kreuz/sourcecode.html. 

Appendix A. Derivation of the SPIKE-distance 

The derivation of the instantaneous SPIKE- measure S{t) 
(Eq. 9) consists of three steps: calculation of instantaneous 
spike time differences, locally weighted averaging, and nor- 
malization. Here we provide additional details about the 
motivation behind each of these steps. 

The first step is the calculation of the absolute instan- 
taneous differences between previous and following spike 
times \Atp{t)\ and |AtF(^)| (Eqs. 4 and 5, respectively). 
By taking into account only the previous and the following 
spike in each spike train the method relies on local informa- 
tion only. The method is also time-scale adaptive, since the 
information used is not contained within a window of fixed 
size but rather within a time frame whose size depends on 
the local rate of each spike train. 

The next two steps, locally weighted average, and nor- 
malization are interchangeable. Here we first calculate the 
locally weighted average of the two differences |Atj(t)|, j = 
P,F: 



(Atj(t)),=p,F = 



E,=P,F|Aij(t)|,/(x|")(i)) 



^j=P,F 



(A.1) 



Using the inverse of the mean intervals from the time in- 
stant under consideration to the previous and the follow- 
ing spikes (xj"-'(t)), j = P,F (averages over Eqs. 6 and 7, 
respectively) as weights we obtain 

\Atp{t)\- 
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which can easily be transformed into 
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(A.3) 

As can be seen in Fig. 1, the sum of the mean intervals 
equals the mean interspike interval: 
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Inserting Eq. A.4 into Eq. A.3 we obtain 

(At,(t)),-=p,F = 



Atp(t)K4"Ht))„ + |AtF(i)|(4"^(t))„ 
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(A.5) 

Finally, in order to achieve time-scale invariance (i.e., all 
stretched and compressed spike trains should yield the same 
value) we divide this locally weighted spike time difference 
by the mean interspike interval. This way we recover Eq. 9: 



S{t) = 



lAtpmxPit))^ + |AtF(t)|(4"^(t))n 



(A.6) 



At the same time the spike time differences are related to 
the local spike rate such that a certain spike time difference 
is the more relevant the higher the rate. This also yields 
to a normalization. The maximum value 1 of S{t) is ap- 
proached when two spikes from the different spike trains 
almost coincide. This can be seen when looking at Fig. IB 
and imagining tp^ getting closer and closer to In this 

case |Atp(i)| approaches a;|g| and |AtF(i)| approaches xjgj. 
Accordingly, S{t) approaches the value 1. Note that overall 
the local increase is outbalanced by surrounding decreases 
so that for converging spikes a reduction of the SPIKE- 
distance is obtained (cf. Section 2.3). Also note that, while 
S{t) can get arbitrarily close to 1 it never actually reaches 
this value, bccaiisc for identical spikes is obtained. Apart 
from the normalization, an analogous derivation holds for 
the multivariate case (Section 2.1.3) as well. For both cases 
the distances are obtained from the instantaneous values 
as their temporal average (Eqs. 10 and 14). 



Appendix B. Previously published measures of 
spike train distance 

Here we restrict ourselves to descriptions of the ISI- 
distance together with its extensions, the Victor-Purpura 
distance, and the correlation coefficient since these are the 
measures that we deal with in more detail. For the other 
measures that we compare against in Appendix D please 
refer to the Appendix of Kreuz et al. (2009). 



B.l. The bivariate ISI-distance and its extensions 

The ISI-distance (Kreuz et al., 2007) and its extensions 
(Kreuz et al., 2009) are based on the instantaneous inter- 
spike intervals. 



B.l. 2. Averaged bivariate ISI-distance 

As in the case of the averaged bivariate SPIKE-distance 
(Section 2.1.2) the averages over pairwise distances and 
over time commute and we can visualize the instantaneous 
average in a time-resolved way. First we calculate the in- 
stantaneous average A{t) over all pairwise absolute ISI- 
ratios \Imnit)\ (cf. Eq. B.l) 
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before we average over time: 



Df = i dtr{t). (B.4) 
Jt=o 

As the original measure, the averaged bivariate ISI-distance 
is restricted to the interval [0, 1]. 



B.l. 3. Multivariate ISI-distance 

We derive a multivariate measure by calculating the in- 
stantaneous coefficient of variation taken across all neurons 
at any given instant in time 
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and averaging over time: 
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For identical spike trains obtains the same lower bound 
of zero as Di and Df, but unlike them it lacks an upper 
bound. 



B.2. Victor-Purpura distance 



B.1.1. Bivariate ISI-distance 

To define a time-resolved, symnicitric. and time-scale 
adaptive measure of the relative firing rate pattern (Kreuz 
et al., 2007) we take the instantaneous ratio between a;|g| 
and a;|gj (Eq. 3), and combine according to: 



h2{t) = 



(t)/xi(i)-l if x[^,{t)<x[^,{t) 



isn 
otherwise. 



(B.l) 

This quantity becomes for identical ISI in the two spike 
trains, and approaches —1 and 1, respectively, if the first 
or the second spike train is much faster than the other. 

The bivariate ISI-distance is obtained by averaging the 
absolute ISI-ratio over time: 

Di = U^ dt\h2{t)\. (B.2) 
^ Jt=o 

The bivariate ISI-distance is bounded in the interval [0, 1]. 



The spike train distance Dy introduced in Victor and 
Purpura (1996) defines the distance between two spike 
trains in terms of the minimum cost of transforming one 
spike train into the other by using just three basic opera- 
tions: spike insertion, spike deletion and spike shift. While 
the cost of insertion or deletion of a spike is set to one, the 
cost cv of moving a spike by some interval is the only pa- 
rameter of the method, and sets the time-scale of the anal- 
ysis. For zero cost, the distance is equal to the difference 
in spike counts, for high costs, the distance approaches the 
number of non-coincident spikes, as it becomes more favor- 
able to delete all non-coincident spikes than to shift them. 
Thus, by increasing the cost, the distance is transformed 
from a rate distance to a timing distance. For multi-neuron 
data, synchrony is defined as the average over all pairs of 
spike trains: 
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B. 3. Correlation coefficient 

In contrast to all the other measures used in this paper 
the bivariate correlation coefficient relies on binning with 

the bin size b (in seconds) as a free parameter. If there is no 
spike within the bin, a zero is assigned, otherwise this bins' 
value is set to one ^ . The two resulting vectors b x and 
b y are used as input for the correlation coefficient C which 
equals the normalized covariance function r according to 

C = C(t^, ty)/^ C(t^, ■ C(ty, 'ty). (B.8) 

Also the multi-neuron correlation coefiicient is calculated 
by averaging over all pairs of spike trains. 

Appendix C. Data 

All recordings and simulations were performed prior to 
and independently from the design of this study. 

C. l. Multi-neuron recordings from retinal ganglion cells 

In Section 3.2.1 the recorded activity of a nearly complete 

population of ON parasol retinal ganglion cells (RGCs) 
is used to compare the performance of the new SPIKE- 
distances with previously published measures such as the 
Victor-Purpura distance and the ISI-distances. The datasct 
used is a reanalyzed subset of the first dataset in Shlens 
et al. (2006, 2009) (cf. Figs. 2A and lA, respectively). 
RGCs were recorded with planar array of 512 extracellu- 
lar microelectrodes, covering an area of 1890 x 900 /xm. 
The spikes of diff'erent cells were isolate (Field et al., 2007) 
and ON parasol cells were identified by their density, their 
response kinetics, and their receptive field characteristics 
(Chichilnisky, 2002). Responses to two kinds of stimula- 
tions were analyzed. Spatiotemporal receptive fields were 
measured using reverse correlation with a white noise stim- 
ulus, composed of a lattice of square pixels updating ran- 
domly and independently over time (Chichilnisky, 2001). 
The intensity of each display primary at each pixel loca- 
tion was chosen independently from a binary distribution 
at each refresh. A stimulus pixel size of 60 /im on a side 
was used. Measurements of spontaneous activity were ob- 
tained in the presence of a spatially uniform, full-field illu- 
mination with intensity equal to the mean intensity of the 
white noise stimulus. For both stimuli, white noise stimu- 
lation and constant illumination, we analyzed 18 segments 
of 10 s each which were distributed regularly over a time 
scan of 30 minutes. 



^ Note that this binary sctuj) of the corrchitioii coefficient disregards 
all differences within the bin size as well as all differences larger 
than the bin size (similar to a rectangular kernel), whereas both 
the Victor-Purpura and the SPIKE-distance weight all spike time 
differences linearly (similar to a triangular kernel). 



C.2. Taste response and temporal coding in the nucleus of 
the solitary tract of the rat 

In Section 3.2.2 we compare the performance of the spike 

train distances in discriminating four different taste stim- 
uli. In Di Lorenzo and Victor (2003) the authors exam- 
ined the reliability of response rate across stimulus repe- 
titions and the potential contribution of temporal coding 
to the discrimination of four taste stimuli (NaCl, sucrose, 
quinine HCl, and HCl) in the nucleus of the solitary tract 
(NTS) of rats. A selected dataset consisting of responses 
recorded from three different neurons is available online 
(http://neurodatabase.org/, Gardner, 2004) as experiment 
"dilorenzo-ndb-1181" . For each of the four stimuli record- 
ings of these cells, labeled #4, #9, and #11, include 19, 23, 
and 16 repetitions, respectively. 

C.3. Hindemarsh-Rose simulations 

The spike trains of the controlled configuration used in 

Appendix D were generated using time series extracted 
from a network of Hindemarsh-Rose neurons (Hindmarsh 
and Rose, 1984). This network was originally designed to 
analyze semantic memory representations using feature- 
based models. Details of the network architecture and 
the implementation of the feature coding can be found in 
Morelli et al. (2005). The clustering properties of these data 
were detailed in Krcuz et al. (2007), their set separation 
was investigated in Kreuz et al. (2009), and we here follow 
the description of these data given in Kreuz et al. (2009). 

In short, the network consisted of 128 Hindemarsh-Rose 
neurons belonging to C/ = 16 different modules with F = 8 
neurons each. The state of the neuron n was determined 
by three dimensionlcss first-order differential equations de- 
scribing the evolution of the membrane potential X„, the 
recovery variable Yn, and a slow adaptation current Z„, 

Xn = Yn-Xl + iXl-Z^ + I^ + an-pn (C.l) 
F„ = 1 - bXl - y„ (C.2) 
Z„ = 0.006[4(X„-1.6)-Z„], (C.3) 

where 

F{U-1) 

Otn= '^nm.Am (C.4) 

m=l 

and 

= yZTi E (C.5) 

m— 1 

are the weighted inter-modular and intra-modular activi- 
ties, respectively, which are dependent on the synaptic con- 
nection weights Wnm and the neuronal activity variables 
An- The external input 7„ was chosen randomly between 
3.0 and 3.1 such that the Hindemarsh-Rose neurons were 
operating in a chaotic regime. 

In a learning stage, input memory patterns were stored 
by updating the synaptic connection weights Wnm between 



different neurons using a Hebbian mechanism based on 
the activity variables An- A neuron was considered active 
whenever its membrane potential X exceeded the thresh- 
old value X — Q and its activity was coded by the variable 
An = 9{Xn — X), where 9 is the Heavyside function with 
9{x) = 1 if a; > and 9{x) = if a; < 0. For these studies we 
restricted ourselves to 29 time series Xn extracted during 
the retrieval stage in which the learned connection weights 
were kept constant. According to their coding properties 
regarding the retrieval of only two distinguished memory 
patterns, they belonged to three clusters: 13 of the neurons 
coded for pattern 1 only, 13 coded for pattern 2 only, and 3 
coded for both pattern ( "Shared" ) . The respective time se- 
ries were labelled by "1" , "2" and "S" followed by an index 
letter. The numerical integration was done using a fourth- 
order Runge-Kutta integration with a fixed step-size of 0.05 
(arbitrary time units). The length of the time series was 
32768 data points. The threshold for spike detection was 
chosen as the arithmetic average over the minimum and 
maximum value of the respective time series. In Appendix 
D.l we follow Kreuz et al. (2007) and evaluate the clus- 
tering of all 29 spike trains, whereas in Appendix D.2 we 
follow Kreuz et al. (2009) and the investigation of the set 
separations was restricted to the two principal clusters. 

Appendix D. Comparison of measures using 
simulated Hindemarsh-Rose time series 

In order to connect with our previous work we finally 
evaluate the performance of the newly proposed SPIKE- 
distances within the bivariate and multivariate setups of 
the simulated network of Hindemarsh-Rose neurons used 
in Kreuz et al. (2007) and Kreuz et al. (2009), respectively. 

D.l. Bivariate setup: Assessing clustering quality 

One important application for measures of spike train 
synchrony is the identification of spike train correlations 
and the recognition of response clusters. In Kreuz et al. 
(2007) we carried out a controlled comparison of six differ- 
ent spike train distances regarding their capability to re- 
produce the clustering within a network of 29 Hindemarsh- 
Rose spike trains (for a description of the data see Appendix 
C.3). 

The measures against which we compared the ISI- 
distance Di were the Victor-Purpura distance Dy (Victor 
and Purpura, 1996), the van Rossum metric Dr (van 
Rossum, 2001), the inversion of the similarity measure 
proposed by Schreiber et al. (2003), and the inverted event 
synchronization Dq (Quian Quiroga et al., 2002). The van 
Rossum metric ZJr measures the Euclidean distance be- 
tween the two spike trains after filtering of the spikes with 
an exponential kernel. The distance Z^sch based on the sim- 
ilarity measure proposed by Schreiber et al. (2003) quan- 
tifies the normalized cross correlation of spike trains after 
Gaussian filtering. The inverted event synchronization 
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Fig. D.l. (color online) Controlled comparison of spike train distances 
on a Hindemarsh-Rose network. A. Comparison of bivariate measures: 
Separation of clusters. B. Comparison of multivariate measures: The set 
separation is defined as the fraction of tests with a Kolmogorov-Smirnov 
statistic equal to 1. The dashed line separates the averaged bivariate 
measures from the multivariate measures. Note that panel A is an 
extension of Fig. 11 from Kreuz et al. (2007), and panel B is an extension 
of Fig. 15 from Kreuz et al. (2009). 

Dq (Quian Quiroga et al., 2002) quantifies the number of 
quasi-simultaneous appearances using a variable time-scale 
that automatically adapts itself to the local spike rates. 
Thus Dq (Quian Quiroga et al., 2002), like the ISI- and 
the SPIKE-distance, does not need a time-scale parame- 
ter, whereas all other distances rely on a parameter that 
sets the time-scale of the analysis. These parameters were 
varied over several orders of magnitude and the highest 
performance obtained was used to represent the measure. 

Details on the analysis can be found in Kreuz et al. 
(2007). In short: We applied the dissimilarity measures to 
all possible pairs of spike trains and from the resulting 
pairwise distance matrices we generated hierarchical clus- 
ter trees. From these dendrograms we extracted the clus- 
ter separation, an indicator which quantifies how well the 
different measures can separate the three clusters of the 
network. We here repeat exactly the same analysis for the 
spike-based distance Ds- As shown in Fig. D.lA both the 
ISI- and the SPIKE-distance almost match the cluster sep- 
aration of the best time-scale optimized measure. The poor- 
est cluster separation is obtained for the optimized -Dsch- 



D.2. Multivariate setup: Assessing set separation 

Following Kreuz et al. (2009) we use the same controlled 
setup with the network of Hindemarsh-Rose neurons to 



compare the performance of several averaged bivariate and 
multivariate measures in distinguishing different levels of 
multi-neuron spike train synchrony. The averaged bivariate 
measures are the extended variants of the bivariate mea- 
sures used in Appendix D.l, multivariate measures com- 
prise the multivariate 1ST and SPIKE-distanccs as well as 
the reliabihties D'^ by Hunter et al. (1998) and by 
Ticsinga (2004). While measures the normalized vari- 
ance of pooled, exponentially convolved spike trains, 
exploits the deviation of pooled spike train statistics from 
the one obtainc^d for a Poisson process. Again the perfor- 
mance of the time-scale dependent measures is optimized 
with respect to the time-scale parameter. 

In order to gradually change the level of synchrony 
in a controlled manner, we construct randomly selected 
sets of spike trains from the two principal clusters of the 
Hindemarsh-Rose network and vary the set imbalance, 
the relative contributions of these clusters. For each set 
imbalance we create five groups of 100 realizations and 
then employ the set separation, a simple measure based on 
Kolmogorov-Smirnov statistics, to quantify how well the 
distributions of measure values for adjacent set imbalances 
can be distinguished. We repeat exactly the same analysis 
for the averaged bivariate SPIKE-distance £>g as well as 
for the multivariate SPIKE-distance f™. Results for all 
measures are shown in Fig. D.IB. 

In this multivariate context the SPIKE-distance Dg 
matches the performance of the best time-scale optimized 
measures. The ISI-distance Df is very close behind. Over- 
all, the averaged bivariate measures are consistently better 
at distinguishing difl[erent set imbalances than the multi- 
variate measures. Among these the highest set separations 
are obtained for the ISI-distance closely followed by 
the SPIKE-distance D™. This is probably due to the fact 
that these measures are not invariant to shuffling spikes 
among the spike trains but rather do take into account the 
spike train of origin for each individual spike. On the other 
hand, measures that act on the pooled spike train, such as 
all measures based on the Peri-Stimulus Time Histogram 
(PSTH), yield the same value regardless of how spikes are 
distributed among the different spike trains (cf. also Fig. 
3). 
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